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We use qualitative arguments combined with numerical simulations to argue that, in the approach 
to the singularity in a vacuum solution of Einstein's equations with T 2 isometry, the evolution at a 
generic point in space is an endless succession of Kasner epochs, punctuated by bounces in which 
either a curvature term or a twist term becomes important in the evolution equations for a brief 
time. Both curvature bounces and twist bounces may be understood within the context of local 
mixmaster dynamics although the latter have never been seen before in spatially inhomogeneous 
cosmological spacetimes. 
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I. INTRODUCTION 



Thirty years ago, Misner [1] and Bclinskii, Khalatnikov and Lifshitz (BKL) [2] noted that Bianchi IX spatially 
homogeneous cosmological solutions of Einstein's equations seem to exhibit a sort of oscillatory behavior in the 
approach to the singularity. This behavior, labeled "mixmaster" [1], involves an infinite sequence of periods (or 
"epochs") during which the solution evolves essentially as a Kasner spacetime [3,4], with each Kasner epoch ended 
by a "bounce" of short duration which changes the evolution from that of one Kasner to that of another one. The 
sequence of Kasners satisfies a rule, called the Kasner map, which takes one Kasner in the sequence to the next. 1 This 
characterization of the Bianchi IX singularity has recently been made rigorous [5] . BKL also made the rather surprising 
claim [2,6] that in spatially inhomogeneous solutions of Einstein's equations, timelike observers 2 approaching a Big 
Bang or Big Crunch singularity should generally see this oscillatory behavior, with the Kasner epoch seen by one 
observer differing from that seen by other neighboring observers, but the sequence of Kasners for each observer still 
satisfying the Kasner map. 

Our previous study of the magnetic Gowdy family of spacetimes [8] provided the first firm support for BKL's claim 
in a spatially inhomogeneous setting. In that work, we numerically evolve spacetimes in the family using the standard 
areal (or "Gowdy") time foliation, and we find that generic observers see oscillatory behavior in the metric evolution. 
Moreover, our studies of the magnetic Gowdy spacetimes indicate that the sequence of Kasners seen by each observer 
follow the pattern of succession predicted by BKL [2,6]. Further, these studies agree with the qualitative picture 
which the Grubisic-Moncrief method of consistent potentials (MCP) suggests [9]. 

Since this magnetic Gowdy work, numerical and MCP studies of two other families of cosmological spacetimes have 
been carried out: the T 2 symmetric vacuum spacetimes and the U{1) symmetric vacuum spacetimes. Both studies 
strongly support the BKL claim that the approach to the singularity is oscillatory. The results for U(l) symmetric 
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1 However, it is not necessarily expected that the evolution converges to a single such sequence of Kasners. It may be that the 
evolution always eventually diverges from any one such sequence and another sequence, which again follows the Kasner map, 
becomes a better approximation. 

2 In this paper we mean by the term "observer" a timelike path with constant spatial coordinates. We assume that a foliation 
and threading have been chosen. Whether results of the sort discussed here will be seen by inequivalent sets of observers is not 
yet generally known. However, this does seem to be true at least in certain cases [7]. 
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solutions have been reported elsewhere [10]. Here, we discuss the behavior near the singularity for T 2 symmetric 
vacuum spacetimes. 

Since U(l) is a subgroup of T 2 = U(l) x U(l), the T 2 symmetric vacuum spacetimes are a subfamily of the U(l) 
symmetric vacuum spacetimes. 3 One may then ask why it is useful to study the T 2 symmetric family directly. The 
reason is that, since the equations for the T 2 symmetric family are considerably simpler (1+1 PDEs rather than 2+1 
PDEs) the numerical studies can be done significantly more accurately. Hence, the studies are more accurate for the 
T 2 symmetric family, and the behavior of the bounces seen by the observers can be monitored more carefully. The 
MCP analysis has been carried out in great detail in this simpler case. We report this study here both to present 
the detailed picture it gives of the dynamics in these spacetimes, and also, we hope, as an aid to obtaining rigorous 
results about the dynamics. 

We define the T 2 symmetric family in Sec. II, noting the relationship between this family and others, such as the 
Gowdy [11] and the Kasner spacetimes. Also in Sec. II, we discuss the areal function and coordinates, recalling results 
which justify their use for T 2 symmetric solutions and writing out the field equations. In Sec. Ill, we set up the MCP 
treatment of the evolution equations for the T 2 symmetric spacetimes and use it to argue that oscillatory behavior 
occurs. We recall that in setting up the MCP form of a given set of evolution equations, one presumes that at each 
spatial point, the fields evolve to Kasner epoch values (not necessarily at the same time for all spatial points); one 
then substitutes these Kasner-like values of the fields into the right hand side of the evolution equations, and attempts 
to infer how the various terms in these equations should behave in time, and what the resulting behavior of the fields 
should be. This analysis predicts that there should be three types of bounces in T 2 symmetric spacetimes: curvature 
bounces, twist bounces and kinetic bounces. A kinetic bounce is not a transition between two distinct Kasner epochs. 
Rather, it occurs within a Kasner epoch. However, in terms of the evolution of the metric functions, it is a bounce on 
a par with the others and its occurrence is necessary for the oscillatory behavior to continue. We state in Appendix A 
the explicit evolution of the fields during each of the three bounces (ignoring in each case terms in the evolution 
equations which are small) and discuss the qualitative nature of each in Sec. III. We compare the MCP predictions 
for bounce behavior with those of BKL. We also discuss in Sec. Ill the MCP argument that, in these spacetimes, an 
observer following a timelike path of constant spatial coordinate should see an unending succession of bounces, a key 
ingredient of mixmaster dynamics and the BKL claims. 

MCP analysis provides useful predictions, but is limited in that, besides being nonrigorous, it does not predict 
whether generic initial data will evolve into a spacetime in which, along each appropriate timelike observer's path, 
a Kasner-like state is reached. (This, again, is a prerequisite for carrying out the MCP study.) To justify the MCP 
predictions, we rely on numerical studies of T 2 symmetric solutions. For representative sets of initial data, Kasner 
epoch values for the fields are reached at each spatial point. Once the Kasner regime is reached at a given spatial point, 
the bounces occur as predicted by our MCP studies, as far as we are able to carry out the evolution. Discussion of these 
results is presented in Sec. IV. These numerical studies do not prove that T 2 symmetric solutions generically exhibit 
oscillatory behavior near the singularity, as predicted by BKL. They do, however, strongly support this contention. 

In the magnetic Gowdy family of spacetimes, we have found that in a generic solution, conditions can occur at 
nongeneric spatial points (e.g., the derivative of a metric component has a zero) with the result that at various points 
near this nongeneric point, there is only a finite number of bounces. While similar conditions occur at nongeneric 
spatial points in a generic T 2 symmetric solution, and while these conditions again appear to persist, there is as yet 
no evidence that the bounces stop at nearby points, in contrast to the situation in a magnetic Gowdy solution. The 
occurrence of the exceptional points is observed in the numerical simulations, and the long time behavior is predicted 
by the MCP analysis. We discuss and exhibit exceptional points in Sections III and IV, but leave extensive discussion 
of these to future work [12]. We make concluding remarks in Sec. V. 

II. THE T 2 SYMMETRIC SPACETIMES 

We define the T 2 symmetric family of spacetimes to consist of globally hyperbolic solutions of the vacuum Einstein 
equations with compact Cauchy surfaces and with a T 2 isometry group acting spatially and without fixed points. 
Generally for spacetimes in this family, at least one of the "twist" functions 

K (X) := fF ,A^rV'I A and K (Y) := e^ pX X»Y»V<>Y x (1) 



3 Note that in studies of the behavior near the singularity in U(l) symmetric spacetimes, a restrictive assumption is made. 
This restriction is consistent with the full range of T 2 symmetric solutions. 
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does not vanish. (Here X and Y are a pair of Killing fields which generate the T 2 isometry group.) If in fact both 
twist functions do vanish, then one obtains the important subfamily of Gowdy spacetimes. 

The Gowdy spacetimes have been extensively studied, and it is believed [9,13-19] that they are all Asymptotically 
Velocity Term Dominated (AVTD). Roughly speaking, this means that as each observer in a given spacetime ap- 
proaches the singularity, she sees at most a finite number of bounces, and eventually settles into a final Kasner epoch 4 
which generally varies from point to point. Since the Gowdy spacetimes are fairly well understood, and since they are 
a set of measure zero in the full family of T 2 symmetric spacetimes, we shall henceforth presume that one or both of 
the twist functions is nonzero, in which case the only topology compatible with T 2 symmetric spacetimes is T 3 x TZ 1 . 

It is very useful in studying the properties of the evolution in a given family of spacetimes to have available a 
universal choice of spacetime foliation which exactly covers the maximal globally hyperbolic development of every 
spacetime in that family. As proven in [21], the "areal foliation" (with corresponding areal coordinates) serves this 
purpose for the T 2 symmetric solutions. We recall that the areal foliation chooses spacelike hypcrsurfaces which are 
invariant under the T 2 action (thereby containing complete orbits of T 2 ) with each leaf of the foliation consisting of 
all orbits of a fixed area. That is, if we let R : T 3 x TZ 1 — > TZ 1 be the function which assigns to a given spacetime point 
the area of the T 2 orbit which contains that point, then the areal foliation chooses for its time function some t cx R. 
In [21] (see also [22]), it is shown that for every T 2 symmetric solution (T 3 x TZ, g) of the vacuum Einstein equations, 
(i) such a function t is indeed timelike, (ii) for every value of t £ (to, oo) with to > (to fixed for each spacetime) 
the t leaf is indeed a T 3 Cauchy surface and (iii) the i-hypersurfaces, with t G (to, oo), collectively cover the maximal 
globally hyperbolic region of (T 3 x 1Z,g). Hence, the areal foliation provides the desired universal choice of time for 
the T 2 symmetric spacetimes. We note that for the Gowdy spacetimes, t is the familiar Gowdy time. 

If we use (x, y) <G T 2 as coordinates labeling points on the T 2 isometry group orbits, and use 9 e S 1 as a coordinate 
parametrizing distinct orbits, then (6,x,y,t) serve as universal coordinates for the T 2 symmetric spacetimes, and we 
may write the generic metric for this family in the form 

g = e^ v - u \-adt 2 + d6 2 ) 

+a e 2U [dx + Ady + (G 1 +A G 2 ) d9 + (M l + A M 2 ) dt} 2 

+a e- 2U t 2 [dy + G 2 d9 + M 2 dt} 2 (2) 

where U, A, v, a, G\, G 2 , M\ and M 2 are functions of 9 and t (independent of x and y), and a is a positive constant. 
This form (2) for the T 2 symmetric metrics is used for the analysis in [21]. Here, to make it easier to compare the 
present study of T 2 symmetric spacetimes with previous similar studies of magnetic Gowdy spacetimes [8] and Gowdy 
spacetimes [13-15,23], it is useful to replace the time function t by t = — hit (one still has an areal type foliation) 
and the metric functions U, A, v and a by the following equivalent functions. 

P = 2U + t, (3) 

Q = A, (4) 

H = — 2 In en , (5) 

\ = Av-AU + 2\na-T. (6) 

In terms of these variables, the metric takes the form 

g = -e^-^^dr 2 + e^ +T ^ 2 d9 2 

+a e p - T [dx + Qdy+(G l +Q G 2 ) d9 + (Mi + Q M 2 ) (-e^dr)} 2 

+cr e- p - T [dy + G 2 d9 + M 2 (-e^dr)} 2 . (7) 

We note, for purposes of comparison, that the metric for magnetic Gowdy spacetimes is the same as (7) except that 
G\ , G 2 , Mi and M 2 vanish. If one relaxes the assumption of the T 2 isometry to allow it to be a local isometry, then 
other spatial topologies in addition to T 3 are possible in the magnetic case or the Gowdy case but not in the general 
T 2 symmetric case with nonvanishing twist. The topology affects the spatial boundary conditions of the functions P 
and Q, but not the qualitative behavior of the evolution toward the singularity [12,16,24]. 

The Einstein vacuum field equations for the T 2 symmetric spacetimes [21] naturally divide themselves into four 
sets. The first set 



4 Asymptotically velocity term dominated behavior is defined more carefully in [16,20]. 
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deK (x) = 



d T K {X) = 



d e K (Y) = 



d T K (Y) = 



(8) 



simply tells us that the twist functions are constant in space and time (and hence are labeled the "twist constants"). 
For any given T 2 symmetric spacetime, we may always replace X and Y by a linear combination of themselves and 
thereby cause one or the other twist functions to vanish (but not both). Hence, without loss of generality, we may 
further presume that only one of the twist constants is nonzero. We label it K. 

The next two sets are the constraint equations and the evolution equations for the metric functions {P, Q, n, A}. 
They involve the twist constant K, but are independent of {G\, G 2 , M\,M 2 }. We discuss these equations below. The 
last set of equations govern {Gi, G 2 , Mi, M 2 }. They take the form 



d T G\ = -e- T d e M x + Q e ^+ 2A + 4P+6 ^/ 4 X, 

d T G 2 = -e- T d 6 M 2 - e (M+2A+4P+6r)/4 Jf- 



(9) 
(10) 



We see from these equations that, once {P, Q, yU, A} have been determined, one obtains {G\, G 2 , Mi,M 2 } by choosing 
Mi(6,t) and M 2 (8,t) to be arbitrary functions of 8 and r, choosing G\(9) and G 2 (9) as arbitrary (initial data) 
functions on S 1 and then integrating Eqs. (9) - (10) over r to get G\{9,t) and G 2 (9,t). Thus {G\, G 2 , M\, M 2 } are 
nondynamical fields. They are essentially "shift functions" , which determine how the coordinates (x, y) evolve in r 
and 9. If K is nonvanishing, {G\, G 2 ,Mi,M 2 } cannot all vanish everywhere in spacetime; the symmetry group does 
not act orthogonally transitively [25] . 

The dynamics of the gravitational field in T 2 symmetric spacetimes lie in {P, Q, fi, A}. To study these fields we find 
it useful to work in Hamiltonian form. Letting Tip, ttq, tt^ and tt\ denote the momenta conjugate to these four fields, 
we find that 7r M may be eliminated, that the functions {P, Q, fi, A, np, ttq, n\} must satisfy the constraint equations 



1 * n 

7T A - -e* = 0, 



(11) 

ITp dgP + TTq 8gQ + 7T\ OgX = (12) 

and that the evolution equations for {P, Q, A, np, ttq, tt\} can be obtained by varying the Hamiltonian density 5 



H 



1 

Att x 



2 , — 2P 2 
TTp + e IT, 



Q 



e- 2T (dgP) 2 + e 2 ^ p - 



\3gQf 



+ aTT X e^ +2P+ ^l 2 K 2 . 



In particular we have 



d T P 

d T TT P 

d T Q 
9 t ttq 



TTp_ 

2tt\ 
1 

2^ 



e- 2P TT 2 Q + e^dggP 



e- 2T (dgP)(d TT X ) 2(P 



(P-t) 



TT\ 



{d e Qf 



-aTT X e (x+2P+ ^ 2 K 2 . 



e- 2P 7T Q 
2tt x 

e 2(P-r) 



2tt x 



[deeQ- {deQ){denx) +2(deP)(deQ)], 



7TA 

■ (2^)2 + e- 2P *Q + e- 2T (dgP) 2 + e 2 ^\dgQ) 2 ] 

+(je (A + 2P + 3r)/2^2 ) 

d T n x = - l -aTT X e^ +2P+ ^l 2 K 2 . 
The evolution for the remaining metric function, ji, follows from the constraint (11) 



(13) 

(14) 

(15) 
(16) 

(17) 

(18) 
(19) 



5 Note that H in (13) is not a superHamiltonian, and is not required to vanish as a consequence of the constraints. It is a 
Hamiltonian (density) corresponding to the choice of time foliation made for these spacetimes. 
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(20) 



The constraint equations, the Hamiltonian and the evolution equations for the fields for magnetic Gowdy spacetimes 
are very similar to these; the main difference is that the twist terms in the Hamiltonian density and in the evolution 
equations are replaced by magnetic terms, with the exponential coefficient for the magnetic terms, e^ x+T ^/ 2 , differing 
from that for the twist terms, e ( A + 2P + 3T )/ 2 . This difference leads to interesting consequences, which we discuss in a 
future work [12]. 

We have already noted the relationship between the familiar Gowdy spacetimes and the T 2 symmetric spacetimes 
discussed here. The (locally) spatially homogeneous subfamily of the Gowdy spacetimes consists of the Kasner 
spacetimes. The (locally) spatially homogeneous subfamily of the generic T 2 symmetric spacetimes, with non- vanishing 
twist, consists of Kasner spacetimes as well. This may seem surprising, since for the standard Kasner Killing vectors, 
all the twist functions vanish. However, one verifies that, in the locally homogeneous subfamily, X and Y (with 
nonvanishing twist) are a linear combination (with constant coefficients) of the three Kasner Killing vector fields. The 
coefficients are constant, but since the norms are changing in time, the angles between the two sets of Killing vectors 
are changing in time. More specifically, consider two orthonormal spatial bases: Ei, made up of eigenvectors of the 
extrinsic curvature (the Kasner directions), and e^, such that each vector in the frame is proportional to a (local) 
Killing vector and such that two of the frame vectors are tangent to the isometry orbits generated by X and Y. Then 
the relation between the two frames is a time dependent rotation. 

Another subfamily of the T 2 symmetric spacetimes is worth noting. If one chooses initial data with Q(6,tq) = 
and ttq(9,tq) = 0, then Q(8,t) = and itq(6,t) = for all points in the spacetime development of this data. (See 
equations (16) and (17).) Hence, one can consider a subfamily, the "polarized" T 2 symmetric spacetimes, with the 
metric coefficient Q — and the corresponding gravitational degree of freedom — turned off. The polarized T 2 symmetric 
spacetimes have been studied using Fuchsian methods, and one finds [26] that there are full-parameter sets of these 
which are AVTD rather than oscillatory near the singularity. Thus although oscillatory behavior is expected to occur 
generically in T 2 symmetric spacetimes, it is not expected to occur in either the Gowdy or the polarized subfamilies. 

III. MCP ARGUMENT FOR OSCILLATORY BEHAVIOR 

The method of consistent potentials (MCP) is a systematic approximation scheme [27] for predicting the behavior 
of cosmological solutions of Einstein's equations in the neighborhood of their singularities. It is based on a key 
assumption, which in practice must be checked numerically. The consequence of this assumption is a weighting of the 
influence of various terms in the Hamiltonian. To describe this, it is useful to split the Hamiltonian density (13) as 
follows: 

H = Hq + Hki n + H sma ii + H curv + Ht W i S t (21) 

where 



H 1 , 2 


(21) 


TJ 1 „-2PJ 

n ~ 4^ Q ' 


(22) 


H small = -Le- 2T (d 9 P) 2 , 
4tt a 


(23) 


H curv = ^-e^ p - T \d e Qf, 
4vr A 


(24) 


H tmst ^a^+ 2P +^' 2 K 2 . 


(25) 




(26) 



The assumption, for a fixed T 2 symmetric vacuum solution (T 3 x 1Z, g) with the singularity at t — > oo, 6 concerns the 
momentary values of the fields for large r. 



6 While the long time existence result [21] does not show that r — > oo, this is expected to be the case generically in this family. 
We need r — > oo to obtain the prediction of an unending sequence of bounces since each (local) Kasner epoch has a finite 
duration in r. 
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Kasner Epoch Assumption (KEA): For each 9 e S 1 , there exists a time T($) such that 

-T(f» « 0, (27) 

-P(0,r (e) )«O, (28) 

P(0,T (e) )-r (e) «0, (29) 

\{9,T (e) ) + 2P{e,T {e) ) + 2,T (e) «0, (30) 

H » fffein, (31) 
i?0 > H smaUl (32) 

-ffo ^ H curv , (33) 
-ffo > Ht W i S t, (34) 

and smc/i t/iai f/ie terms m i/ie evolution equations due to H n dominate the terms in the evolution equations due to 

H-kini H small j Hcurv and H twist- 

Explicitly, this assumption says that for each ^-labeled observer in the spacetime, there is a time T(g) such that all 
of the exponential factors in the Hamiltonian density at (9, T(0)) are very small. (The same factors also appear in the 
evolution equations.) This follows from the conditions (27)-(30). 7 In addition, this assumption says that at this time 
the fields {P, Q, A, 7Tp, ttq, tt\} have developed in such a way that the exponential factors control the relative size of 
the various terms in the Hamiltonian density and also in the evolution equations. 8 

The intent of the KEA is to say that, in any of these spacetimes, the evolution proceeds in such a way that each 
of the 9 = constant observers will, at some time T(g) (generally varying from point to point) reach a Kasner epoch. 
This follows immediately from the evolution equations under the conditions assumed. One might worry that it is 
too restrictive to require that H dominate the other terms in the Hamiltonian density, because it is possible that 
Hq vanishes during a Kasner epoch, but then the prediction of the MCP analysis is that Hq will dominate the other 
terms in the next Kasner epoch occurring at that value of 9, so the KEA will thus be satisfied in this next Kasner 
epoch. 

A consequence of the KEA, combined with the MCP analysis, is that generally (the exceptions are briefly discussed 
later in this Section) the fields evolve in such a way that they do not counteract any explicit exponential decay or 
growth of Hkin, H curv or H tW ist, or that of terms in the evolution equations (14) derived from Hki ni H curv or H tW ist- 
That is, besides dominating the relative size of the various terms in the Hamiltonian density (13) and in the evolution 
equations (14)— (19) , the exponential terms dominate the changes in relative size of the terms. So, for example, if 
\(6,T {e) ) + 2P(6, r ( j)) + 3T (0 ) < and £(A(0,T (fl) ) + 2P(0, t {6) ) + 3T (e) )| (0iT(9)) > 0, then it follows from the KEA that 
Ht W ist is small relative to H at (0, T(0)), and it follows from this assumption combined with the evolution equations 
that H twist is growing in r. 

Numerical results, which we discuss in Sec. IV, indicate that the KEA holds for T 2 symmetric solutions. 

We now consider what happens to the gravitational fields along a fixed 9 observer path after a Kasner epoch has 
begun at some time T($y We presume a fixed spacetime, and we assume the KEA discussed above. Examining the 
evolution equations (14)-(19), we find that the right hand sides of all but (14) and (18) are extremely small. Hence 
the variables {Q, \i, np, ttq, n\} are essentially constant. The variables P and A are not constant; but if we set 

„ :- £ (35) 

we have 

d T P = w (36) 

and 



7 Note that the assumption does not say that the conditions (27)-(30) hold for all r > t^) at 6. 

8 Dominance of the exponential factors is discussed in [8] and assumed in "Assumption A" . There we failed to note that there 
are exceptional situations in which exponential dominance does not hold. We address those cases briefly in this paper and in 
more detail in [12]. 
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d T \ = -w 2 + O 



(37) 



where O indicates terms which, as a consequence of the KEA, can be neglected. The function w is essentially constant, 
so P and A evolve linearly with r. Since all four Hamiltonian potentials, {Hkin, H sma ll, H curv , H tW i S t} are negligible, 
we call the evolution "velocity dominated" when a) the KEA holds, b) {Q, /i, ttp, ttq, tt\} are essentially constant 
(note that this implies that to is essentially constant), and c) 9 r A« —(d T P) 2 = —w 2 . To reiterate, part a) implies 
parts b) and c), that is, the KEA implies that the evolution is velocity dominated at T(gy 

This predicted pattern of evolution for the variables {P, Q, A, fi, ttp, ttq, tt\\ and their spatial derivatives, none 
increasing or decreasing faster than linearly, is consistent with the KEA. The conditions in the assumption continue 
to hold so long as the Hamiltonian potentials stay small relative to Hq. 9 To see whether or not these potentials do 
indeed stay small for r > Trg\ at 9, we need to examine the time derivatives of the exponential quantities in the 
expressions (21) for H kin , H smaU , H curv and H twist . We have 

d T (-2P) = -2to, (38) 
t (-2t) = -2, (39) 
dr{2P-2r) = 2(w- 1), (40) 



Or 



i(A + 2P + 3r) 



1 9 3 

= --to 2 +«,+ -, (41) 



with w approximately constant in r. Clearly the value of the quantity to is crucial in determining whether each of 
the potentials grows or not. In particular, we have, at fixed (6,t^), 

H kin grows if to < (and ttq ^ 0), (42) 

H small decays (if d g P ? 0), (43) 

H curv grows if to > (and dgQ ^ 0), (44) 

Htwist grows if - 1 < to < 3. (45) 

Stating this another way, we have the results listed in Table I. Note that the conditions within the parentheses in 
equations (42)-(44) ensure "generic" behavior. 

It does not follow from Table I that if the KEA holds there must be at least one growing potential at (#, TVm). To 
argue that, we need to assume that both ttq and dgQ are nonzero at (9, T(e))- Since generally these smooth functions, 
ttq and dgQ, are nonzero at a given (6, r), we define (6,t^) to be generic if neither vanishes, and exceptional if 
one of them does vanish. Exceptional behavior (of a different type, to be discussed later) also occurs if to = 1. Our 
subsequent discussion presumes that (6, TYm) is generic unless explicitly stated otherwise; we discuss the exceptional 
cases briefly below and in more detail in [12]. 

Assuming genericity, it does follow from Table I that, during velocity dominated evolution, at least one of the 
potentials is growing. Indeed, the growth is exponential in r. To see what affect this has, we need to consider how 
the fields evolve with one or more of the potentials Hkin, H curv and H tW i S t turned on and hence added to Ho. 

In discussing what happens when some of the potentials become significant, it is useful to keep track of the 
generalized Kasner exponents. These are defined to be the eigenvalues of the extrinsic curvature, divided by the mean 
curvature. It follows from the definition that the sum of the three generalized Kasner exponents equals 1. Defining 
the quantity [9,13-15] 

v:=^-y p + e- 2P TT 2 Q )^ (46) 

the generalized Kasner exponents for a T 2 symmetric solution are 

v 2 - 1 ,„ 2 - 2v , n 2 + 2v m 

tr + 3 v z + 3 v z + 3 

As before, O indicates terms which can be neglected at (9,T(a\) when the KEA is satisfied, as shown in Appendix B. 
Note that the KEA implies v w \ttp\/2tt\ = \w\, and therefore that the generalized Kasner exponents are essentially 



3 Hence the name for this analysis: the "Method of Consistent Potentials.' 
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constant in time. And furthermore, considering the expressions (47), we see that, in addition to K\ + K2 + K3 = 1 
which always is satisfied, when the KEA holds it is also the case that 3 ~ 1. These are the necessary 

and sufficient conditions that a set of three numbers be a set of Kasner exponents. Thus the generalized Kasner 
exponents are approximately a set of Kasner exponents when the KEA holds, and so the KEA does indeed imply that 
the evolution is essentially Kasner at (8,t^). 

We recall that the BKL parameter "u" summarizes the information in a set of Kasner exponents [2,6]. Except for 
the case {0, 0, 1}, there always exists a u > 1 such that 

- ~ M 1 + u "(1 + u) 

Kmin — 1 9 7 Umid — 9 7 Umax — 1 9 l^Oj 

1 + u + u z 1 + u + u z 1 + u + u z 

where K m i n is the smallest of a set of Kasner exponents {K mini n mi ^ n max \, K max is the largest and n m id is the middle 
value. For the expressions (47), the relative magnitude depends on the value of v. for v < 1, one has m < K2 < «3, 
for 1 < v < 3, one has K2 < Ki < K3; and for 3 < v one has Ki < K3 < K\. For these different ranges of values off, we 
obtain different expressions for u: for example, if 3 < v, one has u = (v — l)/2. The expressions for u for each range 
of values of v are given in Table II. 

To understand what happens to the fields when, as a Kasner epoch progresses, one or more of the potentials becomes 
significant, one can look at the evolution of the fields for H = Ho+H k in, H = Ho+Hki n + H tW ist, etc. According to the 
MCP results summarized in Table I, there are five such Hamiltonian densities which in principle must be considered: 
H = H + Hfr in , H = H + H curv , H = H + H twist , H = H Q + H^m + H twist and H = H + H curv + H twist . In 
fact, for both Hk in and H tw i s t, or both H curv and H twist , to become significant simultaneously, some fine tuning is 
needed. For example, if the KEA holds at (80, t^ )) w ith 1 < w < 3, then we have (for constants A » and B ^> 0) 

Hcurv/Ho - e -A+{^-l)(T-T { e a) ) and H twist /H - e -B + (-w 2 +2w+3)(r-T (ea) )/2 _ SQ they &re both growing , TheI1; for 

the given A, B, Ttg \, and the time r' of the next bounce, there is some value of w such that the two potentials are 
equal at r'. A similar argument may be made for — 1 < w < and Hk in /H and H tW ist/H - We find (see Sec. IV) in 
our numerical simulations that simultaneous growth and action of two potentials, though rare, does indeed occur. 



A. Derivation of the Bounce Rules 



If a T 2 symmetric spacetime satisfies the KEA with a particular value of w, and if it approaches one of the five 
types of bounces (H + H kin , H + H curve , H + H tmst , H + H kin + H twist , H + H curv + H tmst ), one would like 
to determine what the MCP predicts for the value of w after the bounce is over. One way to do this is to follow 
the evolution of the fields through each type of bounce (with appropriate Hamiltonian) into the post-bounce Kasner 
epoch, and calculate the change in w directly. The explicit bounce solutions given in Appendix A facilitate this 
approach. Another approach, which we use here, is based on energy and momentum conservation during the bounces. 
That is, using the relevant Hamiltonian density for each type of bounce, and noting that its conservation requires 
certain quantitites to be of the same magnitude — but opposite sign — after the bounce as compared to before, we can 
determine how w changes. Note that in carrying out this approach, it is convenient to consider r(t) to be a dynamical 
variable, dependent on a new time coordinate t. It follows that tt t = —H, and we can work in terms of a Hamiltonian 
constraint density Co rather than H, and treat the former as a superhamiltonian (for unit lapse). For example, for 
H = Ho + Hkin, the Hamiltonian for a kinetic bounce, we have 

Co = = 27r r ttx + X - (4 + 7T 2 Q e- 2P ) . (49) 
For a kinetic bounce, ir T and tt\ are constants of the motion while np changes sign. If we then form the quantity 

— (50) 

we immediately obtain the bounce law for kinetic bounces (where unprimed quantities are evaluated in the Kasner 
epoch before the bounce and primed quantities after the bounce) 

w' = — w. (51) 

We next consider the dynamics for a curvature bounce determined by 

1 



C = = 2tt t tt x + g 



n 2 P + (d e Q) 2 e 2 ^ . (52) 
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Again, ir x is a constant of the motion. If we define 



z = P-r, (53) 

then 

8 t z = tip - 2tt x (54) 
(since the variation of Co yields <9t t = 2ir x ). Now dt z must change sign during the bounce so that 

(ttp - 2n x )' = -(ttp - 2vr A ). (55) 
Dividing both sides by the constant 2tt\ and solving for w' yields the bounce law 

w' = 2 - w. (56) 
We now consider the twist bounce governed by H — H + H tW ist- In this case, we have 

Co = = 2n T tt\ + l -n 2 P + 2a n 2 x k 2 e M+ 2P + 3 ^/ 2 . (57) 

If we evaluate the time derivative of the argument of the exponential in the twist potential using the equations of 
motion obtained from the variation of (57), we find that 

dt{\ + P + y ) = 3n x + n P || + a n x n 2 eS^ 2P+ ^' 2 . (58) 

Asymptotically (i.e. when the twist potential may be neglected), 3ir x + np — TTp/(4:Tr x ) is the momentum associated 
with the time derivative (i.e. growth rate) on the left hand side of (58). Thus, it must change sign during the bounce: 

2 2 

(3n x +7T P - ^L)' = -( 37 r A + 7rp - ^). (59) 

To find a rule for w, we must also recognize (from the equations of motion) that n P — 2ir x is conserved in a twist 
bounce, so that 

(ttp - 2tt x )' = (ttp - 2tt a ). (60) 

If we write both sides of Eqs. (59) and (60) so that a factor of n x and 27r A respectively is shown explicitly, divide 
Eq. (59) by Eq. (60) to cancel the 7Ta's on each side (not the same of course since the left hand side is after and the 
right hand side before the bounce), and identify w, we find that 



2w - w 2 \ /3 + 2w 



w — 1 

Eq. (61) has the solution (there are two solutions but one is trivial) 



(61) 



Note that 1 < w < 3 maps into w' > 3 while — 1 < w < 1 maps into w' < — 1. Thus the former will always yield a 
curvature bounce after the twist bounce while the latter will yield a kinetic bounce after the twist bounce. 

Before deriving the bounce rules for the combined bounces — those with cither H = Ha + Hkin + Htwist or H = 
Ho + H curv + Ht W i S t — we wish to make two observations. We first note that the quantity np — 2n x picks out a 
direction in local minisuperspace which is orthogonal to the twist, so its evolution is essentially unaffected by the 
presence or absence of the twist potential. Secondly, we note that the twist bounce rule (62) may be obtained in a 
way different from that used above. Specifically, we find that the evolution generated by H = H + H tW i S t conserves 
the "energy" 

E = (JJL +an xK 2 e (^+2^+3r)/2 + 3^ ' = ^ +anxK 2 e (A+2P + 3r)/2 + 3^ _ ( g 3) 
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Then, if we factor out ir\ appropriately from Eq. (63), and if we divide on the left and right hand sides by the left 
and right hand sides of Eq. (60), we derive 



II U-l 



(64) 



The nontrivial solution of Eq. (64) yields Eq. (62). 

We now consider the curvature-twist bounce. Analogously to Eq. (63), the evolution generated by H = H + 
H C urv + H twist conserves the energy quantity 

E= lL + (We 2(P - T) +anxK 2 e (A + 2P + 3r)/2 _ np + 5nx _ (65) 
47T A 47T A 

The quantity 7rp — 2ir\ is now not conserved. However, as noted above, the behavior of irp— ir\ during a curvature-twist 
bounce should match its behavior during a curvature bounce; so we have Eq. (55) (as can be verified by considering 
conserved and monotonic quantities). If we now combine Eqs. (55) and (65) as we have described above for Eqs. (60) 
and (63), we derive 



w — 1 / V 1 — w 



w 2 - 2w + 5 



(66) 



This in turn gives the bounce rule 



Note that Eq. (67) maps the region 1 < w < 3 into w < — 1 so that the combined bounce will always be followed by 
a kinetic bounce. 

The previous analysis must be modified for the combined kinetic-twist bounce because successive kinetic and twist 
bounces do not bring w into the range leading to a curvature bounce. To do this, three bounces are required — either 
kinetic- twist- kinetic or twist-kinetic- twist. Either choice leads to the same rule 

W ={—)- (68) 

To avoid any implicit assumption regarding the number of bounces which form the combined bounce, we consider 
only conserved quantities . The energy quantity conserved by H = H + Hkin + Htwist is 

E=^ + + an x e (A+2P+3r)/2 x 2 + 3 (6Q) 

4tt a 4tt a 

The second conserved quantity arises in the event that Hkin = H tW i S t- In that case, irp + 2tt\ is conserved since its 
time derivative is proportional to Hkin — H twist . This leads (using the previous procedures) to 



w 2 + 3\' (w 2 + 3 



w + 1 J \ w + 1 



(70) 



which has the non-trivial solution Eq. (68). Since it is unlikely that Hkin = H tW i S t for any extended time (this would 
require both the growth rates and coefficients to be equal), one would not expect the generic behavior to include 
bounces of this type. In fact, none were seen in the numerical simulations. The bounce rules are summarized in Table 
III. 

Using the bounce laws (51), (56), and (62) and Table II for the BKL parameter u in terms of v = \w\, the change 
in u during each type of bounce may be found. Clearly, Eq. (51) yields u' — u for any kinetic bounce, since the 
change in sign of w does not change v. In a curvature bounce, the initial range w > 1 yields two possible relationships 
between v and u while the final w' can involve all three. The possibilities are shown in Fig. 1(a). All possibilities 
yield the standard BKL map for u: i.e. u' = u — 1 if u > 2 and v! = l/(u — 1) if 1 < u < 2. A similar construction 
for the twist bounce is shown in Fig. 1(b). Here we see that u' = u results for all initial values of w. The rule for u 
is the same in a twist bounce, or in a combination twist-kinetic bounce, as in a kinetic bounce. Furthermore, since 
a solution to any of the three subhamiltonians, H + Hkin, H a + H twist , or H + Hkin + H twist is a one parameter 
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family of Kasner spacetimes, it follows that the generalized Kasncr exponents are (approximately) constant in time 
when any one of these subhamiltonians essentially governs the evolution, and therefore u is (approximately) constant. 
The kinetic and twist potentials are each, at any spatial point, a centrifugal wall. This means [28,29] that the identity 
of the principal axis associated with the growing cosmological scale factor changes during the bounce. However, the 
natures of the bounces are quite different. During a twist bounce or a combined kinetic twist bounce at a point in 
space, the rotation of the principal axes is such that one of them is orthogonal to the T 2 symmetry orbits before 
the bounce and tangent after, while another principal direction is tangent before the bounce and orthogonal to the 
symmetry orbits after. During a kinetic bounce the rotation of the principal directions is in the symmetry plane. One 
principal direction is orthogonal to the symmetry orbit throughout the bounce. In Appendix B a comparison is made 
between the dynamics near the singularity at a spatial point and the dynamics of tilted Bianchi II models studied in 

[30]. _ 

In Sec. IV, we shall examine the validity of KEA and the bounce laws in numerical simulations. However, explicit 
solutions through the kinetic, curvature, and twist bounces are known. These may be used to generate further 
predictions which can in turn be explored in numerical simulations. We leave these for future research. The explicit 
bounce solutions are given in Appendix A. 

B. Exceptional Points 

In a generic T 2 symmetric spacetime there are nongeneric points at which the gravitational field does not evolve 
away from an era of velocity dominated evolution in the manner that we have described thus far. There are three 
cases in which this happens: w = 1 during a Kasner epoch, dgQ = during a Kasner epoch with w > 1, and ttq = 
during a Kasner epoch with w < 0. In each case, one can give rough arguments which indicate that bounces are likely 
to occur. We state these here. 

We first consider the case w = 1. The twist bounce solution given in Appendix A is not defined for w = 1. One can, 
it turns out, write down a solution generated by H = Ho + Htwist explicitly in terms of r in this case. It blows up 
at finite r. It is a Kasner spacetime with dgP = 0, with nonvanishing twist constant, K, and with w = 1. However, 
we claim that this is not a good approximation to the dynamics at a w = 1 exceptional point in a generic spacetime. 
Generally it will not be the case that dgP vanishes. This leads to a situation in which the exponential factors do not 
control the terms they appear in. In the solution that blows up, ir\ — > 0. But if ir\ gets small enough, j^- wrests 

control of H sma u from e~ 2r and H sma u becomes relevant. The subhamiltonian that governs the evolution in this case 
is H = Hq + H srna ii + Ht W i S t- This has (after a canonical transformation) the same structure as the Hamiltonian 
for polarized magnetic Bianchi VIo, in which case there are rigorous results which show that the solution does not 
blow up in finite time, and which predict the bounce rule. One also notes that the t dependence of the argument 
of the exponential in H curv vanishes if w = 1. This means that, generically, H curv contributes a constant (not an 
exponentially decaying) term to Eq. (15) for d T np. But w — 1 yields 7 := np — 2n\ = for the quantity which 
is conserved in the twist bounce. The change in 7 due to the term from H curv will remove the exceptional point 
condition. 

We next consider the case that dgQ crosses zero at some 60 during a Kasner epoch with w > 1. The curvature 
bounce is suppressed in a neighborhood of 6$. The closer to the exceptional point, the longer the bounce is suppressed. 
If w < 2 there will be a twist bounce which sends w — > w' > 3. As the neighborhood on which the curvature bounce 
is suppressed gets smaller and smaller, dggP and ttq grow exponentially if w > 2 [13] (as can be seen in the curvature 
bounce solution (93) by considering the case that £ crosses zero). This wrests control of Hki n and terms in the 
evolution equations derived from H curv from the exponential factors, and causes w to decrease until it is again less 
than two, in which case the twist potential begins to grow again, so bounces will continue. 

In the case that ttq crosses zero during a Kasner epoch with w < a similar mechanism (here dggP and dgQ grow 
exponentially if w < —1) ( [13] and see Appendix A) causes w to become greater than —1, in which case the twist 
potential starts growing, so bounces will continue. 

These arguments suggest that there continue to be oscillations even at or near exceptional points. The behavior of 
the gravitational field at an exceptional point is delicate, however. It may, for example, be the case that higher order 
terms play an important role. Further study of exceptional points is needed, and is now underway. 

C. Minisuperspace picture and twist bounces 

One of the more useful (and most pictorial) ways to study the dynamics of spatially homogeneous cosmological 
solutions of Einstein's equations is via the minisuperspace (MSS) picture. We would like to relate our discussion thus 
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far of local oscillatory behavior and bounces in T 2 symmetric solutions to the MSS approach. In particular, we wish 
to see how the local twist bounces appear in the MSS picture. 

We recall that the MSS approach represents spatially homogeneous spacetimes as follows: For each choice of the 
spatial isometry group G (e.g., G — 1Z 3 <-» Bianchi I, G = SU(2) <-» Bianchi IX) one chooses a fixed group invariant 
frame, and then one can parametrize the set of 3-geometries invariant under G using the MSS variables (volume), 
f3± (anisotropy) and fa (if the metric is not diagonal). Using a simple choice of lapse and shift, one can represent 
a spacetime by a trajectory in the MSS configuration space (Q(t), (3+(t), (3-(t)), with the fa(t) playing a subsidiary 
role. It follows from Einstein's equations, adapted to the spatial homogeneity, that the trajectories corresponding 
to solutions evolve via Hamilton's equations, with the Hamiltonian potential proportional to the scalar curvature 
3 i?(/? + ,/3_,r!). 

For Bianchi I spacetimes, 3 R = 0, so that the trajectories are straight lines. For the diagonal Bianchi IX solutions, 
the potential may be represented by triangular walls, as in Fig. 2. (Note that, in this figure, the dynamics has 
been projected onto the anisotropy plane and rescaled (/3+/|f2|,/3_/|f2|) so that the location of the potential walls is 
independent of fi). The dynamics then consists of straight line segments (effectively Kasner intervals) punctuated 
by intermittent reflections off one or the other of the walls (these are the bounces). The Kasner map describes the 
sequences of bounces, and relates the parameters of the Kasner interval after the bounce to those of the Kasner 
interval before the bounce. The MCP prediction that an infinite sequence of bounces should occur follows from the 
closed nature of the region in MSS configuration space which is bounded by the walls. 

Before tying this picture to the local behavior of T 2 symmetric solutions, we note what happens to the MSS 
description of Bianchi IX spacetimes if the metric is not diagonal. In that case, additional ("centrifugal") walls 
appear, bisecting the angles of the triangle. One or more may be present, depending on the orientation of the 
rotational axis relative to the principal axes in the spacetime. The centrifugal walls do affect the bounces; however, a 
very slightly modified Kasner maps allows one to predict the effect of these bounces on the Kasner interval transitions. 

To obtain a MSS type picture for the local dynamics of T 2 symmetric solutions, it is useful to first do so for certain 
subfamilies. For the Gowdy subfamily (twist equal to zero), we compare the Bianchi I spatial metric 

7/ = e 2i>-40 +d 02 + e 2Q +2/ 3 ++ 2V3/3_ da; 2 + & 2^2^-2^H_ dy 2 (71) 

with the polarized (and diagonal) Gowdy spatial metric 

IPG = e (A+T)/ W + e- T+p dx 2 + e- T - p dy 2 (72) 

to find the following relations between the MSS variables (f2,/3 + ,/3_) and the Gowdy metric components 

P = 2y/S/3- , A = 6(fi-/3+) , r = -2(fi + /3+). (73) 

(These identifications are not unique, since we have singled out one direction, d/dd, in the Bianchi I spacetime to 
identify with the direction of spatial dependence in the Gowdy spacetime.) Then if we rotate in the x-y plane by an 
angle £ and make identifications, t and A are unchanged but we have 

£ P = e 2V3P- CQS 2 £ + £ -2V3P- sin 2 ^ (?4) 

and 

sin £ cos ^e^ 2 ^ 2 ^ 3 - - e 2 "" 4 ^) 

^ ~~ e 20-4/3 + CQS 2 £ + e 2Q+2/3 + +2V3/3- gm 2 £ - < > 

Now adding spatial dependence on 9 to the Kasner solution obtained by the rotation through £ yields a generic Gowdy 
spacetime where we recall that the Gowdy spacetimes may be obtained from T 2 symmetric spacetimes by setting the 
twist constant n to zero and tt\ to 1/2. Unpolarized Gowdy spacetimes have two non-vanishing potentials, H^ in and 
H curv (specialized to the Gowdy case). It follows from (74) that the Gowdy potential H curv will be of order unity 
on the lines labeled C and C in Fig. 2 if that diagram is assumed to depict a local MSS for a Gowdy spacetime. 
Similarly, Hki n will be of order unity on the line labeled K in Fig. 2 when the distances from the system point to C 
and C are equal. 

In a similar fashion, a rotation through an angle lo(t) in the 9-y plane leads to new identifications which correspond 
to those appropriate to a polarized T 2 symmetric spacetime. The identifications are based on the metric (7) with Q 
and Gi set to zero in the spatial metric. We find 

£ -2P = £ -4VS0- CQS 2 w + £ -2V3P- -6/3+ sin 2 ^ (?6) 
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e -2r = e 4/3 ++ 40 ^2 ^ + g 2 ^3/3— 2/3 ++ 4f) ^2 ^ 



(77) 



and the twist potential combination 



(A+2P+3r)/2 _ 



,3/3+ -V3/3- 



2 

- cos w + e 



3/3++V3/3 




(78) 



Each of these identifications contains two exponential terms on the right hand side, only one of which can be large 
at a given time. If we choose one of these to represent the centrifugal wall (e.g. the one consistent with the original 
Kasner identifications for P and r), we find that the centrifugal wall is of order unity along the line labeled T in 



We claim that in the general T 2 symmetric models, the walls C, C , K, T, and the walls allowed by alternate 
identifications (e.g. T" in Fig. 2) are present in the local MSS picture. In that case, the local dynamics is confined 
to the shaded region in Fig. 2. A twist producing rotation yields a term p^/ sinh (6/3+ — 2\/3/?_) in the Kasner or 
mixmaster Hamiltonian [28] . Identification of the exponentials from the metric is not the same as constructing the 
relevant potential — hence, the absence of a minus sign in the denominator of (78). However, if only one of the two 
exponentials in the denominator is large, such a sign cannot be detected. Presumably, an analysis such as that in 
[31] should demonstrate the equivalence between the twist potential and the Kasner centrifugal potential in the local 
MSS picture. 



MCP arguments give a qualitative prediction of what behavior one might see in solutions of Einstein's equations. 
They are, however, neither rigorous nor complete. Thus it is important to compare the solution behavior predicted 
by MCP arguments with the behavior observed in numerical studies of T 2 symmetric spacetimes. As we discuss here, 
the agreement is remarkable. 

Eqs. (14)- (19) are solved numerically using a second order Iterative Cranck-Nicholson (ICN) method (see for 
example [32]). Symplcctic methods used in our previous studies [15,27] fail for these models. Apparently, the operator 
splitting used in the symplectic algorithm allows a pathological behavior in tt\ which is suppressed by other (more 
standard) methods. The CN algorithm can be shown to be more stable but less accurate than symplectic methods 
prior to blow-up of the latter. 

In our numerical studies, we have examined the evolution of the gravitational field for a wide variety of sets of 
initial data. We get qualitatively similar results in all cases. The graphical results displayed and discussed here are 
primarily based on the representative choice of initial data with P = 0, irp = 5 cos(6> + tt/5), Q = cos(6* + 7r/5), ttq = 
0, A = 0, n\ = 1/2 and with a = 1 and k = 1CP 4 . This particular choice of data is useful in that it leads to early 
onset of twist bounces. 

A crucial feature of BKL (or AVTD) behavior is that the evolving fields satisfy the KEA within finite time at 
each spatial point 9. To check for the onset (and later, the recurrence) of field evolution consistent with the KEA, 
we monitor w = TVp/(2n\) at every spatial point. Typical behavior (at three spactial points) is shown in Fig. 3. 
The analysis of Sec. Ill predicts that w(9) should be asymptotically piecewise constant. This is seen to be the case. 
Careful examination (in Fig. 4) shows that the constancy of w(9) becomes an ever better approximation as t — > oo. 

To study the validity of the MCP predictions for the change of the value of w(9) following a bounce, given its value 
before [Eqs. (51), (56), (62), (67), and (68)], we measure the w values in all the Kasner epochs. The next value is 
then predicted using all possible bounce laws. If the new value of w obeys any one of these, then the difference | Atu| 
between the actual and predicted values should be much smaller than those obtained by chance. For the indicated 
initial data, at a typical spatial point, the first twist bounce is followed by a sequence of alternating kinetic and 
curvature bounces. Typically, the comparison of actual to bounce law predictions improves as r increases. Typical 
data are shown in Fig. 5. At some spatial points, a second twist bounce occurs. Fig. 6 shows the values of \Aw\ 
computed for all bounces using the twist bounce law. The small values are indicative of the actual twist bounces. 
Fig. 7 shows the behavior of P at spatial points with a second twist bounce having w > 1 and w < 1. 

Interestingly, we find in our numerical simulations a number of spatial points where a combined curvature-twist 
bounce has occured. This is shown in the plot of combined bounce rule \Aw\s shown in Fig. 8. Fig. 9 shows w(t) 
at such a spatial point with the same quantity at nearby spatial points. At the latter, it is clear that the combined 
bounce has split into separate curvature and twist bounces. 

As illustrated in Fig. 10, for small values of w after the first twist bounce, the evolution changes very little if we 
change the spatial resolution. This indicates convergence (in the computational analysis sense) at such spatial points. 
If w is large after the first twist bounce, then the evolution does appear to be somewhat spatial-resolution dependent 



Fig. 2. 



IV. NUMERICAL RESULTS 
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(as is also illustrated in Fig. 10). This reflects some sensitivity to initial conditions at the first twist bounce where a 
subtraction appears in the denominator of the bounce law. However, within a given trajectory, the agreement with 
bounce laws is more convergent with increasing spatial resolution as seen in Fig. 11. The waveforms at r = 61.66 are 
shown for different spatial resolutions in Fig. 12. Inclusion of adaptive mesh refinement in these codes is in progress 
[33]. 

In Fig. 13, a section of the spatial axis is shown for graphs of dgQ, dgP, and ttq vs at a late value of r. If any of 
these quantities vanish at some 6o, exceptional behavior results there. The zero crossings are shown. For the given 
initial data, the density of exceptional points increases rapidly with time as is shown in Fig. 14. Since such points are 
generated by evolving small scale spatial structure, one could argue that exceptional points should become a dense 
set (of measure zero) as r — > oo. This means that any rigorous statements about the nature of these solutions must 
include consideration of the behavior at exceptional points. 

At a few points, there are anomalies where the bounce laws are violated. This seems to be a consequence of 
inadequate resolution since the effect disappears at higher spatial resolution. The data for P, Q, and A (at two spatial 
resolutions) for the results discussed in this section are shown in Fig. 15. 

The simulations illustrated here cannot be run significantly beyond t w 75 (as indicated on various figures). The 
reason becomes clear from examination of the values of w at r w 75. At several values of 9, w(9) is very close to 
(and less than) unity. From Eq. (62), it is then clear that the next bounce should be a twist bounce with a very 
large new value of w. But w » 1 produces numerical overflows in at least one exponential term (depending on the 
sign of to') in the equations of motion (14)— (19). This is a "physical," resolution-independent numerical instability. 
While possible ways to resolve this problem involve checking the value of w(9) (and thus slowing the code), there is 
no difficulty in principle in either using arbitrary precision arithmetic (see for example [34] ) or the MCP solution (see 
Appendix A) for the next kinetic or curvature bounce. Fig. 16 is constructed by using the twist bounce rule on the 
computed array w(9,t) to obtain w'(9,t). Momentary (pointlike) large values of \w\ arise during bounces when the 
KEA does not hold. The persistent large values of \w'\ for r ~ 75 indicate that dangerously large values are likely to 
arise after the next bounce at some spatial points. 

V. CONCLUSIONS 

We have examined the approach to the singularity in T 2 symmetric vacuum spacetimes. Numerical simulation 
provides strong support for the contention that these models reach an asymptotic regime where the KEA holds. 
Given the KEA, we then predict, using the MCP, the rules relating one Kasner epoch to the next. Again the 
numerical simulations show remarkable agreement with the MCP predictions. 

These spacetimes may be understood as another example family whose members exhibit local mixmaster dynamics 
in the vicinity of the singularity. Yet the local mixmaster dynamics shown here differs from that studied in magnetic 
Gowdy models. In that case, the local MSS potential is closed by a magnetic wall which replaces one of the curvature 
walls that one would expect in a locally Bianchi IX spacetime. In T 2 symmetric spacetimes, the essentially non- 
diagonal centrifugal wall closes the potential. 

Several questions remain open. First, we may ask whether there are an infinite number of bounces. In the absence 
of exceptional points, one could start from any value of w and apply the bounce rules (51)— (62) indefinitely. We 
have argued that the most common exceptional points with dgP = 0, dgQ = 0, or ttq = do not cause the bounces 
to terminate. However, we have seen that the number of exceptional points increases as ever smaller scale spatial 
structure is produced by bounces which occur at different places at different times. Any rigorous discussion of the 
asymptotic behavior of T 2 symmetric models must deal with the exceptional points. We do not yet know the role, if 
any, played by exceptional points where higher derivatives also vanish. Detailed discussions of exceptional points will 
be given elsewhere [12]. 

A second open question concerns the relationship between T 2 symmetric and U(l) symmetric models. If a diagonal 
Bianchi IX metric is expressed in terms of the U(l) variables, all features observed up to now in generic U(l) models 
[10] may be explained in terms of local mixmaster dynamics [31]. It is not yet known whether analogs of the twist 
bounce (a feature of non-diagonal Bianchi IX models) have been missed or suppressed in the existing U (1) simulations. 
A detailed discussion of the relationship between the two classes of spacetimes will be given elsewhere [35] . 

Even with these open questions, we have provided strong support for the validity for T 2 symmetric spacetimes of 
the BKL picture in its most general (local, non-diagonal Bianchi IX) form. We have also provided yet another example 
of the power of the MCP in the analysis of the approach to the singularity in inhomogeneous cosmologies. Finally, 
we have shown how this class of spacetimes allows accurate numerical simulations yet provides a highly non-trivial 
manifestation of local mixmaster dynamics. 
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APPENDIX A EXPLICIT SOLUTIONS FOR THE BOUNCES 



If the gravitational field at (0 O , t^ )) satisfies the KEA, then the local evolution of the gravitational field quantities 
{P, Q, fi, A, 7rp, ttq, 7r\} takes the simple velocity dominated form, but at nonexceptional points at least one of the 
Hamiltonian potentials grows in time. To understand what happens as one of the potentials becomes significant, we 
study the evolution of the fields for each of the three governing Hamiltonian densities H = Ho + Hkin, H = Ho + H curv 
and H = H + H tw i st . Letting {P, Q, fi, A, irp, 7Tq, tt\} denote the data values at (0o,to), with r = T(e ), we obtain 
in each case the explicit 10 general solution on a neighborhood of 0q. 



A.l. Kinetic Bounce 

H = o + Hkin 

= ±-y p + e-^ Q ). (79) 

For a kinetic bounce to occur for some t > To, we need ttq ^ and we need w < 0, or equivalcntly, -kp < 0. We 
presume that both of these conditions hold for the data at (0a, to). 

Now let us define the following series of convenient constants 11 (all depending upon the data at (0 O , r ). 



p:=^-[* P + e-^ Q ]i, (80) 

-P-0T O ~ 

C ■■= \ ^ (81) 

TTp - ZfjlT\ 
e P+0T O 

a ■= TTc^- ^ 

We note that the KEA, with %p < 0, implies that (3 ~ —TTp/2-kx = —w. We also note that the KEA implies that ( is 
very small, and ttq =0<^( = 0. 

In terms of these constants, the solution governed by the density (79) and matching the initial conditions at To is 
(for fixed O ) [13-15] 

P(r) = P - (r - t ) + In ( ^^J , (83) 

Q(r) = Q + <e-^- a{i l ee20Ty (84) 
A(t) = A-/3 2 (t-t ), (85) 



In the case of the twist bounce, we use an implicitly denned function. 
11 These constants depend on 9, and the solution will be a good approximation to the evolution through the bounce on a 
neighborhood of (do)- Knowledge of the solution on a spatial neighborhood is necessary to confirm that the exponential factors 
generically do, both before and after the bounce (and in the case of the neglected potentials also during the bounce), control 
the terms in which they appear. It is also necessary for analysis of the exceptional points, since the spatial derivatives play a 
crucial role. But we are discussing here field evolution in r at the fixed point 9q, so we write the quantities as functions of time 
alone. 
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2 p 2/3r 



1-Ce 
1 + C 2 e 2 ^ T 



7Ta(t) 



7T A . 



(86) 

(87) 
(88) 



A. 2. Curvature Bounce 

H = Ho + Hcurv 

= ^ 2 p + e 2{P - T) (deQn (89) 

In this Hamiltonian, a spatial derivative term appears. However, since ttq does not appear, the general solution 
for fields governed by (89) is relatively straightforward to derive. This Hamiltonian is in fact related to (79) by a 
canonical transformation. The similarity in structure can be seen in the two explicit solutions, but we do not discuss 
this further here. 

For a curvature bounce to occur for some r > r , we need dgQ ^ and we need w < 1, or equivalently, irp < 2n\. 
We presume that both of these conditions hold for the specified data at (9 , t ). 

To write out explicitly the general solution for H in (89), it is again useful to first define a set of constants (depending 
on the data at To): 

P := -5- [(2tt a - itp) 2 + e 2 ^ p -^(deQ) 2 } * , (90) 
2ir\ 

(:=—. ^(1-/3), (91) 

-P+(/5-1)to 

a := = (92) 

1 + C 2 e^ T « 

We note that the KEA, with 2n\ — %p < implies that (i w — (27r A — kp)/2tt\ = (w — 1) and that £ is very small. 
<9eQ = 0^C = 0. The solution takes the following form [36]: 

P (r) = P + (l + /3)(T - To) -l n ^-±^-j, (93) 

Q(r) = Q, (94) 
A( 7 - ) ^A-(l +/3 ) 2 (r-r o) + 21n^ T ±i_j, (95) 

7 r P (r) = 2 7 r A + 2^^--4_j, (96) 

7t q (t) - ttq - e 2 ^"^ {(9 9 a)C - 5(^0 - 2a(9^)Cr + C 2 e 2 <^> [(d e a)( + a(d e Q]} 

e 2 ^ {(d e &)C - a(d e Q - 2a(deP)(r + ( 2 e 2 ^{(d e &)( + a(d e C)}} , (97) 
tta(t) = tt a . (98) 

Note the appearance of P(r) in the exponent of (97); one may substitute in P(r) from (93) to obtain an expression 
in terms of the initial data. 

A. 3. Twist Bounce 

H = H + Ht W i S t 

1 9 \+2P+3t „ 

= 7rp + cr7T A e 5 X 2 . (99) 

47T A 
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For a twist bounce to occur for some t > To, we need —1 < w < 3; that is, —1 < ttp/2tt\ < 3. We presume 
this condition holds. Note that, by assumption, K ^ 0; and note that it follows from the constraint (11) that 
tt\ > 0. Hence H tW i S t is always positive. While we can find the solution to (89) for all values of initial data with 
w e (—1,3), one finds that if w = +1, the solution blows up in finite time [12]. (See Sec. III.) We thus presume that 

we(-i,i)u(i,3). 

As for the other two bounces, in writing out the solution for the twist bounce it is useful to define a few constants 
(depending on the initial data): 12 

S,:=H + H twlst - n P + 5tt\, (100) 

<f> := Htwisu (101) 
v:= K 2 -4(^ P -2^) 2 ]i (102) 

We also find it necessary to work with an implicitly define time: we define the variable T via the implicit equation 

r(T) - t = -1 ln[cosh(i/T) - 8lTX ~ ^ S i n h(^r)] + ^. (103) 

Note that for the case we are excluding, w = 1, t(T) — t = 0. But if w ^ 1, we verify that r(T) is smooth and is a 
strictly increasing function of T, 

Hence, it is invertible, and we can solve for T(t) in principle. (In the numerical calculations which use the symplectic 
algorithm we use a rootfinder to find T(r).) 
If we define 

S(T) = Sinh(i/r) (104) 

K >' v cosh(i/T) - (8tt a - sinh(^T) ' k ; 

the expressions for the solution are the following [36] 

P(t) =P + T-T Q + (n P -2n x )T(T), (105) 

Q(t) = Q, (106) 
A(r) = A — 2 ln[cosh(^T) - (8tt a - smh{vT)/v], 

-4 ln[l - 2<pS(T)] + (2tt a - tt p - 2^)T (107) 

7rp(r) =7T P +4^ A (/)5(r), (108) 

tqW = kq, (109) 

tta(t) = tt a (1 - 2^5(T)). (110) 



APPENDIX B GENERALIZED KASNER EXPONENTS 

When the quantity e = — d T [i/2 is small the generalized Kasncr exponents are approximately 

drA-2e+l 2(u-l) -2(^ + 1) 

9 T A-2e-3' K2 "" 9 r A-2e-3' K3 ~ 9 T A - 2e - 3 ' [ > 

It follows directly from the evolution equations that the denominator appearing in Eqs. (Ill) is bounded above by 
-3. The «'s are exactly the generalized Kasner exponents if the twist constant K = 0. To derive Eqs. (Ill), we first 
define the orthonormal spatial frame, 



The twist bounce solution can be obtained by a canonical transformation of the magnetic bounce solution in the magnetic 
Gowdy case. That solution is given in the context of magnetic Bianchi VIo in [37]. 



17 



m = e -( A +^+-)/ 4 (do -G 1 a x - g 2 d y ), (112) 



e (-P+r)/2 



£/ 2 = ^—d x , (113) 



e (^+r)/2 




U 3 = ^(5 y -Q9 x ). (114) 

The components of the extrinsic curvature in this frame are 

/(0 T A-2e+l)/2 

fc ab = _ e (-^+3-)/4 [ o w_- 1 e p c> T Q ) (115) 

In this frame the twist bounce and kinetic bounce both occur as bounces off centrifugal potentials. It is convenient 
to compare the oscillatory dynamics of the T 2 symmetric spacetimes to that of the tilted Bianchi II models studied 
in [30] using this frame, because the form of the extrinsic curvature is the same in the two cases, and the off-diagonal 
components of (115) are significant during a twist bounce and a kinetic bounce, in turn. The spatially homogeneous 
models studied in [30] have a tilted perfect fluid as source. Note that, without the source, the constraints rule out the 
possibility of oscillatory dynamics in those models, while in the spatially inhomogeneous T 2 symmetric spacetimes 
similar dynamics are obtained in vacuum. Figure II in [30] depicts the oscillatory dynamics which we see at a generic 
spatial point. Figure II i) in that reference depicts the curvature bounce solutions (in our language). Figure II ii) 
depicts the kinetic bounce solutions and Figure II iii) depicts the twist bounce solutions. The identification is fixed 
by setting w — at the point Q 3 in their figure, with w increasing in the clockwise direction around the Kasner circle 
(w — > ±00 at the point T 3 ). The authors of [30] consider an orthonormal frame e , and the variables (E + , £_). Let 
e\ = U2, &i = U3 and e 3 — U\. While the conditions placed in [30] on the spatial frame are not satisfied here, they 
are approximately satisfied at generic spatial points near the singularity. Setting 

_ 1 ( &w-d T \ + 2e-Z \ _ V3 /2™+^A-2e + 3\ 

^ + ~2 \ -9 r A + 2e + 3 )> ~ { -d T \ + 2e + 3 )> {Ub) 

we obtain figure II from [30] by noting that E_/(2 — S+) is constant in a curvature bounce solution, S + — \/3£- is 
constant in a kinetic bounce solution, and S + is constant in a twist bounce solution. 
To continue the derivation of Eqs. (Ill) we next define r = — w and 



a = 



1 if d T Q = and d T P > 

if d T Q = and d T P < ; ( 117 ) 



otherwise 



r V2v 



JO ifd T Q = 0and9 r P = 
b - \ otherwise ' t 118 ^ 



Note that a 2 + b 2 = 1. Consider the orthonormal spatial frame 



The components of the extrinsic curvature in this frame are 



V\ = Ui, (119) 

(-P+t)/2 

V 2 = = — {(a-be p Q)d x + b e p d y }, (120) 

e (P+r)/2 

V 3 = ^{-{aQ + be- p )d x + ad y }. (121) 



, (<9 T A - 2e+ l)/2 u , . „ , , , 
k ab - i e (- A + 3 -)/ 4 ( -b^l v-l ]. (122) 

-a\fl 



2 




If the twist constant K vanishes (so the spacetime is Gowdy) then e vanishes, so the extrinsic curvature is diagonal 
in this frame, and Eqs. (Ill) follow. If the twist constant K does not vanish, the off diagonal components of the 
extrinsic curvature are small except during a twist bounce. 
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The eigenvalues, £j, and eigenvectors, Wi, of the extrinsic curvature are the solutions of k a bW^ = £,ih a bW^ . Per- 
turbation theory for linear operators [38] shows that, when e is small at some point in space, the difference between 
the eigenvalues and the diagonal components and the angles between the eigenvectors and the frame vectors are both 
bounded in terms of e at that point in space. This gives a bound, 50-^/e, for the magnitude of the error in Eqs. 
(111). This bound is not sharp, and holds whether or not the diagonal components of the extrinsic curvature are well 
separated from each other. 

The eigenvectors of the extrinsic curvature are called the Kasner directions, or the principal axes. When e is small, 
the Kasner directions are essentially given by the frame vectors, V{, in that the angle between each frame vector 
and one of the Kasner directions is small. In the solutions of the subhamiltonian Ho + H tW i S t, e grows and decays 
again. We can explicitly compute the rotation of the Kasner directions with respect to the orthonormal frame, Vi in 
the solutions to this subhamiltonian. Note that V\ is orthogonal to the T 2 isometry orbits and the other two frame 
vectors are tangent to the isometry orbits. We find that in each possible twist bounce one of the Kasner directions 
rotates from tangent to orthogonal, and another rotates from orthogonal to tangent. Note that each solution to 
Hq + Ht W i S t is a one parameter family of Kasner spacetimes. In this case we can verify directly that the generalized 
Kasner exponents reduce to a one parameter family of Kasner exponents, constant in time. Since each solution to 
the subhamiltonian H + H^ in + H tW i S t is also a one parameter family of Kasners, it must also be the case that 
during the evolution governed by this subhamiltonian the generalized Kasner exponents reduce to a one parameter 
family of Kasner exponents, constant in time. The bounce rule (71) shows that in this case also, one of the Kasner 
directions rotates from being tangent to the isometry orbits to orthogonal, and a another Kasner direction rotates 
from orthogonal to tangent. 
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FIGURE CAPTIONS 

Fig. 1. Relations between w, w' , u, and u'. (a) Curvature bounce: Initially, w > 1 so that v = w. The dashed 
line shows v' — \w'\ = |2 — w\. Table II is used to compute u for v > 1 and u' for < v'. The horizontal lines show 
1 < u, u 1 < 2. (b) Twist bounce: Initially, — 1 < w < 3. The corresponding u, v, w' , u' , and v' are shown. Note that 
the curves for u and u' = u are superposed. Table II is used to compute u from v. 

Fig. 2. The local MSS in the /3±/|i7| plane. The triangle represents the Bianchi IX MSS potential. Scaling the 
anisotropy variables (3± by the logarithmic volume \ ft\ (fl — ► — oo is the singularity), keeps the bounce locations fixed. 
The Gowdy models begin generically with the walls C and C which disappear as the spacetime becomes AVTD. In 
magnetic Gowdy, a third curvature-like wall M is created by the magnetic field. In T 2 symmetric models, a centrifugal 
wall T closes off the potential. Here the kinetic "wall" is understood to map C onto C so that the dynamics is confined 
to the shaded region. 

Fig. 3. Typical behavior of w(t) at representative values of 9. The upper points are offset by 20 and 40 respectively 
for display convenience. The flat regions characterize the Kasner epochs. 

Fig. 4. Onset and recurrence of the KEA, as monitored by w. (a) w(t,9q) for fixed 9q. Flat regions indicate KEA. 
(b) and (c) Detailed closeups of w(t, 9o) at early (b) and later (c) times. At the later time, w is flatter during the 
Kasner epoch. 

Fig. 5. Typical behavior and accuracy of the bounce laws at three adjacent spatial points. Each point represents the 
smallest difference between a predicted and measured value of w using all of the bounce rules. This shows that each 
sequence starts with a twist bounce and is followed by alternating kinetic and curvature bounces. As the simulation 
evolves, the accuracy of the bounce law prediction improves. The data were obtained by measuring w for all bounces 
over a symmetric region (of length tt) in the simulations. The bounces are numbered consecutively (N) following 
bounces at increasing rata given point and then moving to the sequence of bounces at the next point. The vertical 
lines divide the bounces at a given spatial point from those at the next point. 

Fig. 6. Twist bounces identified by the difference between the measured and predicted values of w. All bounces 
at all spatial points in the considered interval had their preceding and subsequent values of w measured and then 
computed according to the twist bounce rule. Where the difference between the measured and predicted values are 
large, it means that the bounce was not a twist bounce. The twist bounces early in the simulation agree with the 
bounce rule with rather low accuracy because the KEA is not yet completely valid. The highest accuracy agreement 
with predictions indicates second twist bounces later in the simulation. Clustering of the more accurate twist bounces 
just indicates that similar behavior is occurring at nearby spatial points. 

Fig. 7. Behavior of P at twist bounces. Since in a Kasner epoch, d T P ~ w, the piecewise constant w implies a 
piecewise linear P. The twist bounces are indicated by the arrows. If Wo > 1, then w' > 3 and the next bounce will 
be a curvature bounce. If wo < 1, then w' < — 1 and the next bounce is kinetic. 

Fig. 8. Combined curvature and twist bounces. This graph is generated as in Fig. 6 but using the rule for combined 
bounces. The actual combined bounces have |Au>| m 10~ 6 . 

Fig. 9. Structure of a combined bounce. The combined bounce's w(t) at 9o is shown as a solid line. The segment 
after the bounce has w < —1 so a kinetic bounce will follow; w(t) for a point with 9 slightly less than 9o is shown 
with the dotted line. Here the final prc-kinctic bounce segment (K) is preceded by pre-twist (T) and pre-curvature 
(C) segments. The dot-dashed line shows w(t) for 9 slightly greater than #o- Here first a pre-curvature bounce and 
then a pre-twist bounce segment precedes the final pre-kinetic bounce segment. 
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Fig. 10. Spatial resolution dependence. The evolution (for w(6, r)) is shown at three representative values of 9 
(offset by 15 and 30 respectively) for 1024 (solid line) and 2048 (dotted line with circles) spatial grid points. The 
dependence on spatial resolution increases with the value of \w\ which follows the initial twist bounce. 

Fig. 11. | Aw | for a sequence of bounces at the same value of 9 for different spatial resolutions. The difference 
between the measured and predicted values of w is shown. 

Fig. 12. Resolution dependence of waveforms. As has been noted in the evolution of U(l) symmetric cosmologies 
[10], narrowing spiky features [15] cause the simulations to yield resolution dependent results where the functions are 
not smooth. The choice of initial data made here yields an especially spiky waveform for P. A representative portion 
is shown. 

Fig. 13. Exceptional points. Exceptional points with deQ = and dgP = are associated with the peaks in P 
while ttq — causes apparent discontinuities in Q [15,13]. Zero crossings of all three functions are shown at a late r 
value for a portion of the #-axis. 

Fig. 14. The number of exceptional points vs r. The growth in the number of exceptional points vs r is shown. 
While N appears to level off, this could just reflect the exponential increase of Kasner epoch duration characteristic 
of mixmaster dynamics. N(t) is not a power law. 

Fig. 15. P(9,t) (top), Q(9,t) (middle), and X(9,t) (bottom) are shown for the full simulation (with arbitrary 
scales for their values). The left hand column uses 1024 and the right 2048 spatial grid points. In each frame, the 
horizontal axis is — n/5 < 9 < 9n/5 and the vertical axis < t < 76. 

Fig. 16. Limits on the simulation. The plot shows w'(9,t) computed from the simulations w(9,t) using the twist 
bounce rule. The scale is set so that values > 100 ( < —100) appear white (black). The white and black lines which 
extend to the end of the simulation indicate ^-values which are destined to have dangerously large values of w after 
the next twist bounce. Only a portion of the 9 axis is shown. 
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TABLE I. Behavior of Potentials when KEA holds. 



Condition 


Potentials that Grow 


Potentials that do not Grow 


w < -1 


and 




H kin 




H small 


H C urv Hfrwist 


-1 < w < 


and 




ttkin H-twist 




H small 


Hcurv 


< w < 1 


H-twist 


Hj^in 


H small 


Hcurv 


1 < w < 3 


and 




Hcurv Htwist 


H kin 


H small 




3 < iy 


and 


<%Q/0 


Hcurv 


H kin 


H small 


Htwist 



TABLE II. Calculation of BKL parameter u. 



Range of v 


< v < 1 


1< w < 3 


3 < v 


Value of u 


u= i±» 

1 — If 


« = A- 

if— l 





TABLE III. Summary of bounce rules. 



Bounce type 


Kinetic 


Curvature 


Twist 


Curvature- Twist 


Kinetic- Twist 


Bounce rule 


w' = — w 


w' = 2 - w 
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FIG. 1. Relations between w, w' , u, and u . (a) Curvature bounce: Initially, w > 1 so that v = w. The dashed line shows 
v' — \w'\ — \2 — w\. Table II is used to compute u for v > 1 and v! for < v'. The horizontal lines show 1 < u, u' < 2. (b) 
Twist bounce: Initially, — 1 < w < 3. The corresponding u, v, w' , v! , and v' are shown. Note that the curves for u and v! = u 
are superposed. Table II is used to compute u from v. 
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FIG. 2. Fig. 2. The local MSS in the (5±/\Sl\ plane. The triangle represents the Bianchi IX MSS potential. Scaling the 
anisotropy variables /3± by the logarithmic volume |fi| (SI — > — oo is the singularity), keeps the bounce locations fixed. The 
Gowdy models begin generically with the walls C and C' which disappear as the spacetime becomes AVTD. In magnetic Gowdy, 
a third curvature-like wall M is created by the magnetic field. In T 2 symmetric models, a centrifugal wall T closes off the 
potential. Here the kinetic "wall" is understood to map C onto C so that the dynamics is confined to the shaded region. 
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FIG. 3. Fig. 3. Typical behavior of w(t) at representative values of 9. The upper points are offset by 20 and 40 respectively 
for display convenience. The flat regions characterize the Kasner epochs. 
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FIG. 4. Onset and recurrence of the KEA, as monitored by w. (a) w(t, 8q) for fixed 9o. Flat regions indicate KEA. (b) and 
(c) Detailed closeups of w(t, 6q) at early (b) and later (c) times. At the later time, w is flatter during the Kasner epoch. 
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FIG. 5. Typical behavior and accuracy of the bounce laws at three adjacent spatial points. Each point represents the smallest 
difference between a predicted and measured value of w using all of the bounce rules. This shows that each sequence starts 
with a twist bounce and is followed by alternating kinetic and curvature bounces. As the simulation evolves, the accuracy of 
the bounce law prediction improves. The data were obtained by measuring w for all bounces over a symmetric region (of length 
71") in the simulations. The bounces are numbered consecutively (TV) following bounces at increasing r at a given point and 
then moving to the sequence of bounces at the next point. The vertical lines divide the bounces at a given spatial point from 
those at the next point. 
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FIG. 6. Twist bounces identified by the difference between the measured and predicted values of w. All bounces at all spatial 
points in the considered interval had their preceding and subsequent values of w measured and then computed according to the 
twist bounce rule. Where the difference between the measured and predicted values are large, it means that the bounce was 
not a twist bounce. The twist bounces early in the simulation agree with the bounce rule with rather low accuracy because the 
KEA is not yet completely valid. The highest accuracy agreement with predictions indicates second twist bounces later in the 
simulation. Clustering of the more accurate twist bounces just indicates that similar behavior is occurring at nearby spatial 
points. 




FIG. 7. Behavior of P at twist bounces. Since in a Kasner epoch, d T P — w, the piecewise constant w implies a piecewise 
linear P. The twist bounces are indicated by the arrows. If wq > 1, then w' > 3 and the next bounce will be a curvature 
bounce. If wo < 1, then w' < —1 and the next bounce is kinetic. 
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FIG. 9. Structure of a combined bounce. The combined bounce's w(r) at 8o is shown as a solid line. The segment after the 
bounce has w < — 1 so a kinetic bounce will follow; w(t) for a point with 6 slightly less than 8o is shown with the dotted line. 
Here the final pre-kinetic bounce segment (K) is preceded by pre-twist (T) and pre-curvature (C) segments. The dot-dashed 
line shows w(t) for 9 slightly greater than 8q. Here first a pre-curvature bounce and then a pre-twist bounce segment precedes 
the final pre-kinetic bounce segment. 
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FIG. 10. Spatial resolution dependence. The evolution (for w(8,t)) is shown at three representative values of 6 (offset by 
15 and 30 respectively) for 1024 (solid line) and 2048 (dotted line with circles) spatial grid points. The dependence on spatial 
resolution increases with the value of Itul which follows the initial twist bounce. 
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FIG. 11. \Aw\ for a sequence of bounces at the same value of 6 for different spatial resolutions. The difference between the 
measured and predicted values of w is shown. 




FIG. 12. Resolution dependence of waveforms. As has been noted in the evolution of U(l) symmetric cosmologies [10], 
narrowing spiky features [15] cause the simulations to yield resolution dependent results where the functions are not smooth. 
The choice of initial data made here yields an especially spiky waveform for P. A representative portion is shown. 
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FIG. 13. Exceptional points. Exceptional points with dgQ = and dgP = are associated with the peaks in P while ttq = 
causes apparent discontinuities in Q [15,13]. Zero crossings of all three functions are shown at a late r value for a portion of 
the #-axis. 
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FIG. 14. The number of exceptional points vs r. The growth in the number of exceptional points vs r is shown. While 
N appears to level off, this could just reflect the exponential increase of Kasner epoch duration characteristic of mixmaster 
dynamics. N(t) is not a power law. 




FIG. 15. P(8,t) (top), Q(6,t) (middle), and X(9,t) (bottom) are shown for the full simulation (with arbitrary scales for 
their values). The left hand column uses 1024 and the right 2048 spatial grid points. In each frame, the horizontal axis is 
— 7r/5 < < 9n/5 and the vertical axis < r < 76. 
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FIG. 16. Limits on the simulation. The plot shows w'(8,t) computed from the simulations w(8,r) using the twist bounce 
rule. The scale is set so that values > 100 ( < —100) appear white (black). The white and black lines which extend to the end 
of the simulation indicate ^-values which are destined to have dangerously large values of w after the next twist bounce. Only 
a portion of the 9 axis is shown. 
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